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Abstract The authors present a new simple algorithm to approximate 
weakly stochastic differential equations in the spirit of |10||17|. They ap- 
ply it to the problem of pricing Asian options under the Heston stochastic 
£> \ volatility model, and compare it with other known methods. It is shown 

that the combination of the suggested algorithm and quasi-Monte Carlo 
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1 Introduction 

1.1 The Problem and its Motivation 

We consider a stochastic differential equation written in the Stratonovich 
form 

r-i d pt 

Y(f,x) = x + Vo(Y(s,x)) ds + Y Vi(Y(s,x))odB[, 

Jo ti J « (1) 

VjeC; (r n ;R n ), 

where B = (b\- ■ ■ ,B d ) is a standard Brownian motion, and C™ (r n ;R n ) 

denotes the set of R N -valued smooth functions defined over R N whose 
derivatives of any order are bounded. In particular, we will use the clas- 
sical notation Vf(x) = Z%i V 1 (x) (df/dxj) (x) for V e q°(R N ; R N ) and / a 
differentiable function from R" into R. This stochastic differential equation 
can be written in ltd form: 

Y(t, x) = x + V ( Y(s, x)) ds + Y \ V i ( Y ( s ' x )) dB i> 
Jo j~f Jo 

where 

1 d 

7=1 

Now, given a function / with some regularity, how can one approxi- 
mate efficiently E [f ( Y(l, x))] ? It is equivalent to the following deterministic 
problem: if L is the differential operator V + (1/2) £ i=1 Vj and u is the so- 
lution of the heat equation 

— (t, x) = Lu, u (0, x) = f(x), 

how does one approximate u (1, x) (which is equal to E [f ( Y(l, x))] by 
Feynman-Kac theorem |8|). 

This problem has had a lot of attention because of its practical impor- 
tance: it gives the evolution of the temperature in some media, and also 
represents price of financial derivatives under stochastic financial models 
such as Black-Scholes 1 1 1. 

Non-probabilistic methods to solve the PDE (such as finite difference 
methods) seem to only work well when L is elliptic and in low dimension. 
We refer to [14J for a more detailed discussion on the subject. We will focus 
in this paper on probabilistic methods. 
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1.2 Notation 

If V is a smooth vector field, i.e. an element of C£°(r n ;R n ), exp(V)x 
denotes the solution at time 1 of the ordinary differential equation 

^7T = V fe) / z = x. 

For x e R, |_xj denotes the integer part of x. For a random variable X, Var[X] 
denotes the variance of X. 



1.3 Probabilistic Methods 
1.3.2 Order 1 

The most popular probabilistic method to approximate E [/ (Y(l, x))] is 
called the Euler-Maruyama method |9|. We first fix n independent d- 
dimensional random variables Z\, ■ ■ ■ , Z„ such that, if X denotes a standard 
normal random variables, 

E[p(Z,)] = E[p(X)] (2) 

for all polynomial p of degree less than or equal to 3. Then one defines 
recursively the following random variables: 

Y (EM),« _ 

,A.q X f 



(EM) „ = (EM),„ 1 ? / (EM),„X + J_fy. ( X (EM),nX * 



Then, one can show 1 9 1 1 25 1 that for an arbitrary C 4 function / 

||E[/(xf M) '»)]-E[/(Y(l /a :))]||<C / i. (3) 

Of course, one needs an algorithm to compute E [/(Xj EM ''")j . If the Zj are 
constructed from Bernoulli random variables, E^/^X^ BM ^'"^j is a discrete 

sum, but one would need to do 2 nd additions, which can be rather lengthy 
when nd is large (one is then forced to do some Monte-Carlo on a discrete 
measure). If the Z& are normal random variables, one then is forced to do 
use some Monte Carlo or quasi-Monte Carlo techniques. When nd is big, 
quasi-Monte Carlo method become less effective than Monte-Carlo, but if 
nd is not too high, quasi-Monte Carlo method can be very efficient. 

Another method with the same rate of convergence appeared in 1171 . 
and is called cubature on Wiener space of degree 3. It is defined with the 
following recursive formula: 



Y (cub3),n _ 
A Q — X, 



f - d 



^-(cub3),n 
k/n 
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Such algorithm can be seen as a practical application of the Wong-Zakai 
theorem 1 8 1 1 26 [, when the Zfc are normal random variables. 
If B" = (B"' 1 , . . . B"' d ) (n e N) is the piecewise linear approximation of the 
Brownian motion defined by 

B'/ = flnfj + 1 - nf)B LnfJ/n + (nt - \nt]) B (L „ t j +1) /„, 

and Y" denotes the solution of the ordinary differential equation 

Y» = x + f Vo (Y?) & + Y" f Vi (Y») dBf 1 , 
Jo " Jo 

then the Wong-Zakai theorem states that Y" converges almost surely to 
Y T . It is easy to see that x^ cub3 ^'" and Y" are equal in law, proving the 
convergence of the weak algorithm cubature on Wiener space of degree 3 
(but this argument does not provide the rate of convergence). 

Remark 1 In the algorithm cubature on Wiener space of degree 3, one has to 
solve numerically ODEs (unless one is lucky and one has a close form solu- 
tion!). One possibility is to take its Taylor approximation of order 1 for the 
approximation of exp ( V) x and we fall back on an Euler scheme. Taking a 
better approximation (Taylor approximation of order 2) will give a scheme 
sometimes described as the Milstein scheme. Not spending enough care 
on the approximating method of the ODEs to be solved can result in some 
catastrophic situations. A general case where that happens is when the 
diffusion is almost surely on a subset of K", that is, does not fill the whole 
space. If one has an approximation scheme which at some time provides an 
answer outside this set (which is what happen if one approximates badly 
the ODEs), the algorithm may go very wrong or even bug. Increasing n 
(which is costly) or artificial techniques can be implemented to solve this 
problem, while this can be overcome by taking an appropriately good ap- 
proximation of the ODEs which have to be solved (we usually recommend 
a high order Runge-Kutta scheme, or an adaptive step size scheme, but 
this may depend on the particular SDE to approximate). We will give an 
example of this problem in Section[3] 

Remark 2 Random variables which satisfy are easy to find. One can take, 

for a fixed i, Z. to be d independent Bernoulli or Gaussian random variables. 
A more elaborate choice of such random variables appeared in [17J[22J. 

Remark^ Here, we have used the subdivision (k/n) k€ \ 0r . n } of [0,1] . It is 
not clear whether taking equal time steps is optimal or not. Recently, 
Kusuoka 1 12 1 proved that the partitioning into equal time steps is opti- 
mal when we use the algorithm which we will propose in this paper. We 
do not want to address this problem in this paper, and we will always take 
subdivisions with equal time steps. 
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2.3.2 Higher order 

A way to obtain approximations of higher order is based on the under- 
standing of more terms in the stochastic Taylor formula (see 1 3 1 and 1 9 1 for 
example). When the vector fields V, commute, it is relatively easy to find 
a scheme of high order, see 1 9 1 and the references within. In the general 
case, one needs to understand how to approximate weakly the increments 
of the Brownian motion together with its first few iterated integrals. This 
was first successfully done, to our knowledge, in |10||16||23||24||13| and 
then generalized with the method cubature on Wiener space 1171 . 



1.4 Romberg Extrapolation 



Consider a nice scheme of order p, that is, a scheme Xj!™ d ^'" such that for 
smooth /, there exists a constant Kf such that 



E[f{x^ n )]-E[f(Y(l,x))]-K f ^ 



< C 



/„p+r 



Then, 



2P 



2P 



i4/( x ! ordp, 1]-^i E [/( x !° [dp l 



(4) 



provides a scheme of order p + 1. We refer once again to |25| for more 
details and the proof that the Euler-Maruyama scheme and its successive 
Romberg extrapolations are "nice" schemes. Recently, it was proved that 
our new algorithm presented below is a "nice" scheme 1 12l . 



1.5 A remark on the Monte Carlo method 

Let W be a random variable. When we compute E[W] by Monte-Carlo 
method with M samples, we consider a random variable (Tj^Li Wkj /M 
where W,-'s are independent random variables whose distributions are 
identical to Ws. We denote this random variable by MC( W,M). By virtue 
of the central limit theorem, we can consider that MC (W,M) behaves as a 
normal random variable of mean E[W] and variance Var[ W]/M. 

Let xj° r p ''" denotes a scheme of order p of the type above. To calculate 
^(,ordp),n numer j ca ]iy^ one neec i t approximate an integral over a nC(d) 
dimensional space (C(d) denoting a function depending on d; for Euler or 
Cub3, C(d) = d. As we will see later, C(d) = d + 1 for our new algorithm). 
If one uses the Monte-Carlo method to approximate this integral, and uses 

M samples, the random variable MC(/(Xj° rd, ^'"),M) is considered. The 
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situation is summarized as following relations: 

E [/(Y(l, x))] =E[f (Xf rdp) '")] + O (nn , (5) 

Varf/(x; ordp) '")]^ 



MC (f(x { ° rdp) '"),M) ~ N E f/(xf rd?,) '")l , 



M 



(6) 



Two types of approximation errors are involved in this calculation. One 
is the difference between E [f(Y(l,x))] and E [/(Xj 0rdp ''")J and the other is 

the difference between MC (/ (xf ldp) ' n ) , m) (co) and E [/ (x< ordp) '")]. In this 
paper, we call the former error discretization error and the latter error inte- 
gration error. shows that we can consider the integration error of Monte 
Carlo method to be a normal random variable of mean and variance 
Var[/(x< ordp) '")]/M. 

Because the difference between Var [/(x< ordp) '")] and Var [/(Y(l,x))] is 
very small, we get the following remark. 

Remark 4 As long as we use the Monte Carlo method for numerical ap- 
proximation of E[/(Y(1, x))], the number of sample points needed to attain 
the given accuracy is independent of the dimension of integration, namely 
the number n of partitions and the order p of the approximation scheme. 



1.6 A remark on the quasi-Monte Carlo method 

Although there are some results which justify the quasi-Monte Carlo 
method and give theoretical error with respect to the number M of sam- 
ple points and the dimension of the integral domain, those results help 
little for error estimation in practice when we apply the quasi-Monte Carlo 
method to weak approximation of SDEs (see |20| or |21[). The following 
observation seems to be widely accepted: 

Remark 5 In contrast to the Monte Carlo case, the number of sample points 
needed by the quasi-Monte Carlo method for numerical approximation of 
E[/(Y(l,x))] depends heavily on the dimension of integration. Smaller the 
dimension, smaller number of samples are needed. 

The integral that we have to approximate to obtain x| ordp ''" is on a space 
of dimension nC(d). If the numerical method is of high order and nC (d) 
is not too big, one can then use quasi-Monte Carlo with this numerical 
method to obtain a very fast algorithm. 

Therefore, it seems optimal to look for a (simple) scheme of order greater 
than that of the Euler-Maruyama scheme (one), with C (d) remaining com- 
parable to d (i.e. the C (d) of the Euler-Maruyama scheme). This is the object 
of this paper, where we suggest a new numerical scheme of order 2, with 
C (d) = d + 1. We will show its efficiency by numerically pricing an Asian 
option under the Heston model. 
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2 Presentation of the new Algorithm 

We present our new algorithm, of order 2. 

Theorem 1 Let (A,Zi) ie{lr .^ be n independent random variables, where each 
Ai is a Bernoulli random variable independent of Z;, which is a standard d- 

dimensional normal random variable. Define |X^ ew) '"} i t =0/ ... < „ to be a family of 

random variables as follows: 

-^(New),n _ x 
A (fc+l)/n ~ 



exp g) exp 






f 4 v A 


:f: 


•••exp 












exp g) exp 




•••exp 





exp 



exp 



/ Vq\ y(New),« 

\2n) k >» 



(New),n 
k/n 



if A* = +1, 
if A = -1. 



(7) 



Then, for all f £ C^°(1R N ), 

\E[f(x^")]-E[f(Y(l,x))] 



that is, our new algorithm is of order 2. 
A few remarks before all: To compute 



exp 



exp 



<Z\V X 



yjn t 



■ exp 



<z{v A 



exp 



X 



(New),n 
k/n 



one needs to solve d + 2 ordinary differential equations. First along the 
vector field V from t = to t = 1 1 {In) with starting point X[^ ew) '*, then 

along Vd from t = to £ = Zt / yfn with starting point the solution of the 
ODE we have just solved, and we repeat similar operations d + 2 times. 
One would need an algorithm to solve this ODE numerically (unless one 
has a close form solution), and we, once again, strongly suggest that one 
pays a lot of attention to the quality of such algorithm. 

One of course will have to use an algorithm to approximate E [ / ^ New) '" ) j , 
but this is just a (difficult but classical, common to Euler algorithm for ex- 
ample) problem of integrating a function on a finite dimensional space. The 
simplest but quite effective method is to do some basic Monte-Carlo sim- 
ulation of the random variables (A,-, Z,) 



ie{l,- ,n 



. One could also simulate the 



random variables (A;, Z;), e n ... n \ with some quasi-Monte Carlo techniques, 
or replace the random variables Z, with some discrete random variables 
with the right moment up to order 5. As this is a very classical problem and 
common to all the other probabilistic solutions to our numerical problem, 
we do not provide anymore precisions here. 
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Proof The proof is quite classical, so we will not go into details. The reader 
should be convinced that the algorithm is of order 2 once we show that for 
/ smooth enough, 



\E[f(^")]- E umii-,x))] 



n 5 



The error over n steps, from the Markov property of Y, would then be n 
times n~ 3 . We consider a smooth function /. First observe that, from the 
Feynman-Kac theorem, 



E [/ (Y(l/n, x))] - (/(x) + hf(x) + ^E 2 / (x)) 
Developing E 2 , that means 



< C' f n- 3 . 



/(X)+ I L/(X) + _L L 2 /(X)=/(X) + I 



d \ 

' v} 



1=1 



/(x) 



In 2 



i=i 



/(x). 



Now we need to approximate E [/ (x^ w) '")] . Using Taylor approxima- 
tion of the ODEs involved, we quickly see that the absolute value of 

/(exp(^Vo)ex P (^y 1 )...exp(^y d )exp(-Ly ), 



minus 



/(*) + 

1 

2^2 



1 d 



Vf 



+ 



1=1 

d 



/(X) 



i=l i=l i=l i<; 



/(X) 



is bounded by C"n 3 . Inverting the order in which the vector fields are 
integrated, we obtain that the absolute value of 

/(exp(iy )exp(iz^ (( )...exp(iz ; v,)exp(i ; y ) J 



minus 



f(x) + 



1 d 



+ 



2n 2 



d A 

1=1 
d 



/(*) 



vfvj 



1=1 



1=1 



i=l 



/(x) 
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is bounded by C"n 3 . Adding up and dividing by 2, we obtain that 



|E[/(x^ ew) '")]-E[/(Y(l/n^))] 



< 



C + C" 



Remark 6 Using the results in 1 10 1 and [11J, one can show the convergence 
of the algorithm with / Lipschitz continuous, under a condition on the 
vector fields weaker than Hormander condition. We do not do it here to 
avoid writing a very technical paper. 

This algorithm could be seen in a non-trivial way as a particular case 
of the algorithm cubature on Wiener space of degree 5. One should also 
notice some common features with splitting methods. 



3 Numerical Example: Application to Finance 

In this section, we numerically compare our new algorithm to the Euler- 
Maruyama scheme and their Romberg extrapolation. We calculate the price 
of an Asian call option with maturity X and strike K written on an asset 
whose price process Y\ satisfies the following two factor stochastic volatil- 
ity model (Heston model \7\): 



Y 1 {t,x) = x 1 + f \iYi{s,x)ds+ f Yi(s,x) ^Y 2 (s,x)dB 1 (s), 
Jo Jo 

Y 2 (t,x) = x 2 + f a{0- Y 2 (s,x)) ds + f S ^Y 2 (s,x)dB 2 (s), 
Jo Jo 



(8) 



where x = {x\,x 2 ) e (K >0 ) 2 , (B 1 (f), B 2 (f)) is a 2-dimensional standard Brown- 
ian motion, and a, 6, ji are some positive coefficients such that 2a - S 2 > 
to ensure the existence and uniqueness of a solution to our SDE |5[. The 
payoff of this option is max (Y^T, x)/T - K, 0), where 

Y 3 (t,x)= f Yi(s,x)ds. (9) 
Jo 

The price of this option becomes DxE [max (Ys(T, x)/T — K, 0)] where D is 
the appropriate discount factor. We set T = 1, K = 1.05, jU = 0.05, a = 2.0, 
= 0.1, 6 = 0.09, and {x lr x 2 ) = (1.0, 0.09). We ignore D in this experiment. 

Let Y{t,x) = \Yi{t,x),Y 2 {t,x),Y 3 {t,x)). We transform the SDEs (SJ and © 
into a Stratonovich form SDE: 



2 r t 



Y(t,x)= > , Vi(Y(s,x))odB'(s), (10) 
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where 

(' (yi, 1/2, 3/3)) = I yi \fi - y ) , a(6 - y 2 ) - —, yi 

Vi( f (yi,y2,y 3 )) = '(yi^o / o) ( n ) 
v 2 ('(yi,y2,y3)) = '(o,^v^'0). 

3.1 Implementation of the algorithm 

We apply the algorithm which we introduced in Section|2]to this problem. 
3.1.1 Solutions of the ODEs 

We can easily get exp (sVi) and exp (sVy (s € R) as follows: 

exp (sVi) f (y X/ y 2 , y 3 ) = (yie s ^ y 2 , y 3 ) , 
exp(sy 2 )'(yi,y2,y 3 ) = yiJy+Vyz) ,y& • 



(12) 



As there exists no closed form solution to exp (sy )/ we are forced to use 
an approximation and we choose: 



exp (sV Q ) f (yi, y 2 , y 3 ) = W S X ^(s), g 3 (s)) , (13) 



where 



gl (s) = y l exp^- J ~) S + ^(e- as -l) 
g2(s) = J + (y 2 -J)e-" s , 

g3( S ) = y 3 + _L_ Uo(s 3 ), 

J=6-^-, and A = u-^. 

The error compared to the true solution is O (t 3 ^ in small time t, creating an 

additional error of O (w~ 3 ) at every step of the algorithm, but as the error of 

our scheme at every step was also O (ft -3 ) / taking the above approximation 
of exp (sVq) does not alter the convergence rate of the algorithm. 

Following the same discussion, it is easy to see that we have to approx- 
imate exp(fVo) i n such a way that the order of the produced erro r is 0(£ 4 ) 
when we use Romberg extrapolation, which we introduced in ll.4l together 
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with the algorithm. In this experiment, we approximate gs(s) by the tradi- 
tional order 4 Runge-Kutta method when we use Romberg extrapolation. 

Here, we see that one of the advantages of this algorithm over the 
Euler-Maruyama scheme is the one we mentioned in Remark^ When we 
apply the Euler-Maruyama scheme to this process JS}, it may happen that 

the square volatility process ( Y 2 )[ EM ^'' 1 becomes negative, and the algorithm 
then fails at the next step (as we will have to take its square root). On the 
other hand, equations $12\ and fL4|l show that our new algorithm does not 
share this problem. There exists a way of avoiding this problem with the 
Euler-Maruyama scheme HI . 

3.2.2 A remark on general implementation 

In general, it is not always possible to obtain the closed form solution 
to exp(sV^). Even in such cases, it is not difficult to implement our new 
algorithm. All we have to do is to find an approximation of exp(sy ) whose 
error is 0(s 3 ) and approximations of exp(s V,), (z + 0) whose errors are 0(s 6 ). 
This can be achieved by Runge-Kutta like methods and we can find some 
examples of them in |2- 

We remark that when we use Romberg extrapolation together, we have 
to approximate exp(sV ) with 0(s 4 ) error and exp(sy ; ) (z + 0) with 0(s 7 ) 
error. 

3.2.3 Application of the quasi-Monte Carlo method 

Our new algorithm has the virtue that the application of the quasi-Monte 
Carlo method to this algorithm is possible in a straight forward way, once 
we embed (A,-, Z,-)^., .,„) into [0, \) n{ - d+l \ This is an advantage of the algo- 
rithm over algorithms proposed in 1 18 1, 1 19 1, and 1 13 1 which also enable us 
to proceed higher order weak approximation. 

3.2 Comparison to Euler-Maruyama scheme 

We compare numerically our new algorithm to the Euler-Maruyama scheme 
with and without Romberg extrapolation. Such methods involve, as we 
saw, approximation of an integral over a finite dimensional space; we will 
do these approximations using the Monte Carlo method and the quasi- 
Monte Carlo method. 

There are many studies on acceleration of Monte Carlo methods 1 6 1 but 
we choose the crude Euler-Maruyama scheme with and without Romberg 
extrapolation as only competitors by the following reasons: 

1. Only our new algorithm and the Euler-Maruyama scheme are very 
universal and applicable easily to any type of problems described in 
subsection ll~T1 

2. Almost all of variance reduction techniques which we can apply to the 
Euler-Maruyama scheme are also applicable to our new algorithm. 
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Discretization Error and Num. of Partitions 



0.01 



0.001 



1e-04 



1e-05 



1e-06 



1e-07 



Euler-Maruyama 
Euler-Maruyama + Extrpltn 
New 

New+ Extrpltn 
0(1/rtf 
0(1/n J ) 
5e-05 
5e-06 



10 



100 1000 
Num. of Partitions 



10000 



100000 



Fig. 1 Error coming from the discretization 



These are important advantages of our new algorithm. Many existing al- 
gorithms lack one or both of these properties. For example, in 1 15 [, they 
proposed the trapezoidal algorithm which accelerates Monte Carlo pricing 
of Asian option price. But this algorithm works only for pricing of Asian 
option written on one dimensional diffusion. There are many such type of 
problem-specific algorithms and we exclude them, because in this paper 
we focus on universal algorithms which work for any weak approximation 
problem of any diffusion processes defined by Q. 
In this experiment, we consider 

E [max (Y 3 (T,x)/T - K, 0)] = 6.0473907415 X 10" 2 

which is obtained by our new algorithm with extrapolation, quasi-Monte 
Carlo, n = 96 + 48, and M = 8.0 x 10 9 . 



3.2.2 Discretization Error 

Figure ^ shows the relation between the number of partitions in our dis- 
cretization of the interval [0, 1] (n in the description of the algorithm) and 
the error of the algorithms. We observe that to achieve 10~ 4 accuracy, our 
new method with Romberg extrapolation requires n = 6, our new method 
needs n = 12, while the Euler-Maruyama scheme with Romberg extrapola- 
tion needs n = 24, and the simple Euler-Maruyama scheme needs n > 2000. 
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In all algorithms, consumed time is proportional to n X M, where M is the 
number of sample points. 

3.2.2 Convergence Error from Monte Carlo 



QMC : New + Extrpltn, 1/At=4 — i- 
QMC : New, 1/At=12 — x- 
QMC :Euler-Maruyama + Extrpltn, 1/At=16 — *■ 
MC : Euler-Maruyama s- 

5e-05 

5e-06 




1e-06 




1000 10000 100000 
Num. of Sample points 



1e+06 



1e+07 



Fig. 2 Convergence Error from quasi-Monte Carlo and Monte Carlo 



We have already mentioned in II. 51 that the convergence performance of 
the Monte Carlo method is independent of the number of partitions. We 
can see in Figure|2|that in this experiment this statement holds. This figure 
also shows that to achieve 10~ 4 accuracy with 95% confidence level (2a) by 
using Monte Carlo method, we need over 10 8 sample points. We can also 
see in this figure that the Monte Carlo errors which come from algorithms 
boosted by the Romberg extrapolation become greater than those of the 
original algorithms. 

3.2.3 Convergence Error from quasi-Monte Carlo and Monte Carlo 

Figure |3 also shows that the performance of the convergence of the quasi- 
Monte Carlo method depends on the number n of partitions and on the 
algorithms. Figure |2] seems to show that the quasi-Monte Carlo method 
outperforms the Monte Carlo method specially when used with our new 
algorithm and that the algorithm needs 2 X 10 5 sample points for 10~ 4 
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accuracy, the algorithm with extrapolation 2 X 10 5 sample points, and Euler- 
Maruyama with extrapolation 5 X 10 6 sample points when we use the 
quasi-Monte Carlo method. 

3.2.4 Performance comparison with respect to consumed time 



Table 1 #Partition, #Sample, and CPU time required for 10 4 accuracy. 



Method 


#Partition 


#Sample 


CPU time (sec) 


E-M + MC 


2000 


10" 


1.72 x 10 b 


E-M + Extrpltn + MC 


16 + 8 


10* 


2.06 x 10 j 


New + MC 


12 


10* 


1.24 x 10 3 


New + Extrpltn + MC 


4 + 2 


10 s 


6.2 x 10 2 


E-M + Extrpltn + QMC 


16 + 8 


5x10" 


1.28 x 10 2 


New + QMC 


12 


2x 10 b 


3.3 


New + Extrpltn + QMC 


4 + 2 


2x 10 b 


1.73 



The elapsed time of each method required for 10~ 4 accuracy is shown in 
Tabled We find in this table that our new algorithm with Romberg extrap- 
olation and the quasi-Monte Carlo method provides the fastest calculation. 
Our new algorithm with Romberg extrapolation and quasi-Monte Carlo is 
about 80 times faster than Euler-Maruyama scheme with Romberg extrap- 
olation and quasi-Monte Carlo. We also see that even without Romberg 
extrapolation, our new algorithm is still faster than any boosted Euler- 
Maruyama method. 

At last we would like to mention Remark |4] and Remark |5] again. The 
remarkable performance of our new algorithm is closely related to the 
property of the quasi-Monte Carlo method noted in Remark|5J 
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